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(Dated: February 8, 2008) 

The theory of adiabatic invariants has a long history and important applications in physics but 
is rarely rigorous. Here we treat exactly the general time-dependent 1-D harmonic oscillator, q + 
uj 2 (t)q = which cannot be solved in general. We follow the time-evolution of an initial ensemble 
of phase points with sharply defined energy Eq and calculate rigorously the distribution of energy 
E± after time T, and all its moments, especially its average value E\ and variance fi 2 . Using our 
exact WKB-theory to all orders we get the exact result for the leading asymptotic behaviour of fi 2 . 

PACS numbers: 05.45.-a, 45.20.-d, 45.30.+S, 47.52. +j 



Adiabatic invariants, usually denoted by /, in time 
dependent dynamical systems (not necessarily Hamilto- 
nian), are approximately conserved during a slow process 
of changing system parameters over a long typical time 
scale T. This statement is asymptotic in the sense that 
the conservation is exact in the limit T — > oo, whilst for 
finite T we see the deviation AI = If — Ii of final value 
of If from its initial value li and would like to calculate 
AI. Here we just remind that for the one-dimensional 
harmonic oscillator it is known since Ehrenfest that the 
adiabatic invariant for T — oo is I = E /lo, which is the 
ratio of the total energy E = E(t) and the frequency 
of the oscillator u>(t), both being a function of time. Of 
course, 2~kI is exactly the area in the phase plane {q,p) 
enclosed by the energy contour of constant E. A general 
introductory account can be found in [lj and references 
therein, especially 0,0- However, in the literature this 
AI is not even precisely defined. As a consequence of 
that there is considerable confusion about its meaning. 
Let us just mention the case of periodic parametric res- 
onance in one-dimensional harmonic oscillator where the 
driving is periodic and yet e.g. the total energy of the sys- 
tem can grow indefinitely for certain system parameter 
values. (In this work we give a precise meaning to these 
and similar statements.) Therefore to be on rigorous side 
we must carefully define what we mean by AI. This can 
be done by considering an ensemble of initial conditions 
at time t — just before the adiabatic process starts. Of 
course, there is a vast freedom in choosing such ensem- 
bles. In an integrable conservative Hamiltonian system 
the most natural and the most important choice is tak- 
ing as the initial ensemble all phase points uniformly dis- 
tributed on the initial iV-torus, uniform w.r.t. the angle 
variables. Such an ensemble has a sharply defined initial 
energy Eq. Then we let the system evolve in time, not 
necessarily slowly, and calculate the probability distribu- 
tion P(Ei) of the final energy E\, or of other dynamical 
quantities. 

This is in general a difficult problem, but in this work 
we confine ourself to the one-dimensional general time- 
dependent harmonic oscillator, so N = 1, described by 



the Newton equation 



q + LU 2 (t)q = 



(1) 



and work out rigorously P{E\). Given the general depen- 
dence of the oscillator's frequency uj(t) on time t the cal- 
culation of q(t) is already a very difficult unsolvable prob- 
lem. In the sense of mathematical physics is exactly 
equivalent to the one-dimensional stationary Schrodinger 
equation: the coordinate q appears instead of the proba- 
bility amplitude if), time t appears instead of the coordi- 
nate x and u; 2 (t) plays the role of E — V(x) = energy - 
potential. In this paper we solve the above stated prob- 
lem for the general one-dimensional harmonic oscillator, 
but the details of our calculations are delegated to an- 
other publication £|. 

We begin by defining the system by giving its Hamilton 
function H = H(q,p,t), whose numerical value E(t) at 
time t is precisely the total energy of the system at time 
t, and for the one-dimensional harmonic oscillator this is 



(2) 



where q,p,M,LO are the coordinate, the momentum, the 
mass and the frequency of the linear oscillator, respec- 
tively. The dynamics is linear in q,p, as described by 
(|T|). but nonlinear as a function of tu(t) and therefore is 
subject to the nonlinear dynamical analysis. By using 
the index and 1 we denote the initial (i = to) and final 
{t = t\) value of the variables, and by T = t\ — to we de- 
note the time interval of changing the parametrs of the 
system. 

We consider the phase flow map (we shall call it tran- 
sition map) 



5o 
Po 



(3) 



Because equations of motion are linear in q and p, and 
since the system is Hamiltonian, $ is a linear area pre- 
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serving map, that is, 



a b 
c d 



(4) 



with det($) =ad — bc= 1. Let E'q = H(q 0} p , t = t ) be 
the initial energy and E\ = H(qi,pi,t = ti) be the final 
energy, that is, 



Ei = l( iCq0+ J P ° )2 + M.l(a qo+bP of). (5) 

Introducing the new coordinates, namely the action / = 
E/ui and the angle 0, 



Expression is positive definite by definition and 
this leads to the first interesting conclusion: In full gen- 
erality (no restrictions on the function u>(t)\) we have 
always E± > Eqlui/ujq and therefore the final value of 
the adiabatic invariant I± = E\/u)x is always greater or 
equal to the initial value Iq = Eq/luq. In other words, the 
value of the adiabatic invariant never decreases, which is 
a kind of irreversibility statement. Moreover, it is con- 
stant only for infinitely slow processes T = oo, which is 
an ideal adiabatic process, i.e. /i = 0. For periodic pro- 
cesses oj\ = ujq w e see that always E\ > Eq, so the mean 
energy never decreases. The other extreme to T — oo is 
the instantaneous (T = 0) jump where ljq switches to oj\ 
discontinuously, whilst q and p remain continuous, and 
this results in a = d = 1 and b = c = 0, and then we find 




cos 0, Po = V2ME a sin (6) 



from (J3J) we obtain 

Ei — E (acos 2 + 0sin 2 + 7 sin 20), (7) 

where 



2 ? 



+ a 2 -%, = d l +u;tM z b z , 

M*Wq UJq 

7=t^-+^M— • (8) 

Given the uniform probability distribution of initial an- 
gles equal to l/(2-7r), which defines our initial ensemble 
at time t = 0, we can now calculate the averages. Thus 



Ei = -^I E 1 dcf>=^-(a + 0). 
That yields E\ — E\ = E (5 cos 20 + 7 sin 20) and 



El 



^ = i Ei-Eif = ^-(S 2 + 1 2 



(9) 



(10) 



where we have denoted S = (a — 0)/2. 

It follows from (JSJ, |JS} that we can write (|TT)j) also in 
the form 



» 2 = (Ei-Eir = ^ 



Ei\ 2 w? 
w 



S 



(11) 



It is straightforward to show that for arbitrary positive 
integer m, we have (Ei — Ei) 211 ^ 1 = and 



(Ei -Eif 



(2m- 1)!! 



({Ei - Ex)' 



(12) 



Thus 2m-th moment of P(Ei) is equal to (2m — 
l)!!/i 2m /?7i!, and therefore, indeed, all moments of P(Ei) 
are uniquely determined by the first moment Ei . 



F - Eo (^ 



Ei 



1)> M = "IT 



^-1 

2 
^0 



(13) 



Below we shall treat the special case with cj 2 = 2cl> 2 , and 
thus will find fi 2 /E 2 = 1/8 = 0.125. 

Our general study now focuses on the calculation of the 
transition map namely its matrix elements a,b,c,d. 
Starting from the Hamilton function J2J) and its Newton 
equation Q we consider two linearly independent solu- 
tions ipi(t) an( i fait) an d introduce the matrix 



* {t> \Mipi (t) MM*) 



(14) 



Consider a solution q(t) of such that 

«(*o) = «>, h(h)=Po/M. (15) 

Because "0i and ^2 are linearly independent, we can look 
for g(i) in the form 

g(t) = A^i(t) + B^ 2 (t). (16) 
Then A and 5 are determined by 

U)=*M:)- (i7) 

Let qi = q(t x ), p x = M<j(*i). Then from fT5 ft — (fTT ^ we 

see that 



91 ] =*(i 1 )*- 1 (to)( 90 
Pi J \Po 



(18) 



We recognize the matrix on the right-hand side of 118|) 
as the transition map <I>, that is, 



<T> -\ a c b d )=V(ti)y-\t ). 



(19) 



Due to lack of space we mention only one specific model, 
namely the linear model defined by the piecewise linear 
co 2 (t) function 



io 2 (t) = { uj 2 + M^I) t if 0<t<T 



if t < 

if < t- 
if t > T 



(20) 



3 



In this case the equation can be solved exactly in 
terms of the Airy functions, and the formalism explained 
above leads in a straightforward but lengthy manner to 
the final exact result for E\ and consequently for /i 2 etc. 
It is too complex to be shown here. The special case 
uJq = 1 and lo\ = 2 has been checked very carefully, also 
numerically, and y? goes correctly from 1/8 at T — to 
zero as T — * oo, in a typical oscillatory way. Using the 
well known asymptotic expressions for the Airy functions 
we find the leading asymptotic approximation 



E 2 



(e 1 - E 1 y 

E 2 




4\/2 cos(- 



8 a/2, 



3e 



(21) 

where we introduce the adiabatic parameter e = 1/T 
which is assumed small. Please observe the oscillatory 
approach to zero as e — * 0, which in the mean goes to 
zero quadratically as e 2 . 

Returning to the general case we now mention that the 
final energy distribution function written down as 



i-V 

2vr^ 



do 



dE 1 



(22) 



<t>=<i>j (Ei) 



cannot be calculated analytically in a closed form in any 
useful way, because it boils down to finding the roots of 
a quartic polynomial, so we do not try to do that here, 
although numerically it shows interesting aspects. It has 
a finite interval as its support, between the lower limit 
E m i n and the upper limit E max , and at both values it has 
an integrable singularity of the type 1/^fx. In between 
for every value of E\ = const — E\ (fa , this equation has 
four solutions, namely fa, fa, fa, fa, and thus we have 
to sum up all four contributions in the general formula 
(|22|l . On the other hand, as we have seen, the moments 
of this interesting distribution function can be calculated 
exactly to all orders. 

We proceed with the calculation of the transition map 
$ in the general case, and because Q is generally not 
solvable, we have ultimately to resort to some approxima- 
tions. Since the adiabatic limit e — > is the asymptotic 
regime that we would like to understand, the applica- 
tion of the rigorous WKB theory (up to all orders) is 
most convenient, and usually it turns out that the leading 
asymptotic terms are well described by just the leading 
WKB terms. 

We introduce re-scaled and dimensionless time A 

A = et, e = 1/T (23) 
so that JIJ is transformed to the equation 



Let q+(X) and q_ (A) be two linearly independent solu- 
tions of (|24l) . Then the matrix (|14H takes the form 



( q+ (A) «_(A) 

\eMq' + (X) eMq'_{\) 



(25) 



and taking into account that Ao = eio, Ai = tt\ we obtain 
for the matrix 1|19[) the expression 



<I> r \\ =* A (A 1 )^ 1 (Ao). 



(26) 



We now use the WKB method in order to obtain the 
coefficients a, b, c, d of the matrix $. To do so, we look 
for solution of (|24|l in the form 



q(X) = wexp < —cr(X) 



1 



(27) 



where cr(A) is a complex function that satisfies the differ- 
ential equation 



(ct'(A)) 2 W(A) = -u 2 (X) 



(28) 



and w is some constant with dimension of length. The 
WKB expansion for the phase is 



a(X) = J2^kW- 



(29) 



fc=0 



Substituting (|29|l into (|28|l and comparing like powers of 
e gives the recursion relation 

-. n— 1 

a' 2 = -c 2 (A), a' n = -ttt(Y. °i<-k + <-i)- (3°) 



fe=i 



Here we apply our WKB notation and formalism [j| and 
we can choose a' + (A) = iw(A) or <j' _(A) = — iw(A). 
That results in two linearly independent solutions of l|24|) 
given by the WKB expansions with the coefficients 

cto,±(A) =±i / u(x)dx, <7!,±(A) = -jrlog 

i /- A 3^(x) 2 -2u{x)u"{x) 
(7 2 , ± =±- / —3 dx, ... (31) 

Since u>(X) is a real function we deduce from H3Q(I 
that all functions cr' 2k+1 are real and all functions o~' 2k 
are pure imaginary and a' 2k , = — o' 2k °2k+l + = 
a 2k+i - wnere k — 0, 1, 2, ... , and thus we have a' + = 
A(X) '+ ifl(A), a'_ = A(X) - iB(X) where A(X) = 

Er=o^ +1 <+i(A), B(X) = -iE^o^k+W. 
Integration of the above equations yields 



V'(A) + ^ 2 (AMA)=0. 



(24) 



<j+ =r(X) +is(X), a- = r(A) -is(A), (32) 
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where r(A) = J Xg A(x) dx, s(A) = £?(x) efe. Below 
we shall denote si = s(Ai). 

Using this notation we find that the elements of the 
transition matrix $ have the following form, after taking 
into account that det($) = ab — cd = 1, 



Hence, in the case n = 2k — 1 we can assume 

A(X) =e 2k - 1 cr' 2k _ h+ (X) + h.o.t. 
B(X) =uj(X) - ie 2k a' 2k + {\) + h.o.t. 



(36) 



1 



b = 



V B Q Bi 
1 



Aq sin ( —) — B cos ( — 



sin — 

M^BqBI V e 



A I 



(A A 1 



IhUh i siu 



V BqBi 

{A B t - A 1 B )cos(^) 

A lS in(^) +S lC os(^) 



(33) 



vBqB 



This is so far exact result, based on the WKB expan- 
sion technique. What we are mostly interested in is the 
asymptotic behaviour of /i 2 when e is small and tends 
to zero. All other aspects and technical details will be 
published in a separate paper Q. 

Let us consider the first order WKB approximation, 
that is, 



A(A) « eoi )+ (A), B(\)t 
We find for the variance 



(A) 



w(A). (34) 



M 2 2 


/ 2/2 


/2 








/I 4 C0S 






0(e 3 ). (35) 



and obtain 



4fc _ 2 / ^2fc-l.+ (Al)- ^2fc-l.+ ( A o) 

W l (T 2fc-l,+ ( A o) cr 2fe-l.+ ( A l) /2si 
1 5 : COS 

^0 V e 

+ (9(e 4 ' £ - 1 ). (37) 
In the case when n — 2k we can suppose 



A(X) =e 2 «+ l a 2k+lt+ (\)+ h.o.t 
B(X) =u(\)-ie 2k a 2k)+ (\)+h 

Then, similarly as above, we obtain 



(38) 



E 2 



o- 2 fc,+ ( A i) . ^2fc, + ( A o) 



2i0 2 



From this we can conclude that if w(t) is of class C m 
(having m continuous derivatives, m — n — 1) /z 2 goes 
to zero oscillating but in the mean as oc e 2n = g 2(m+i) 
If m — oo (analytic functions) according to Landau and 
Lifshitz the decay to zero is oscillating and on the 
average is exponential oc exp(— const /e). 



Substituting into the last formula w(A) = vT+A we ob- 
tain exactly the approximation lj21|l . 

Suppose now that all derivatives at Ao and Ai van- 
ish up to order (n — 1), i.e. oj'(Xq) — uj'(Xi) = • • • = 
u( n -V(\ ) = w( n - x )(Ai) = 0, and lo^ {X )lu^ (Ai) ^ 0. 
Then ai(Ao) = ^(Ai) = •■■ = ^(Ao) = <_i(Ai) - 
0, (t;(A K(A 1 ) ^0. 
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